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BACKGROUND OF THE INVENTION 

1. Field of the Invention 

[0003] The present invention relates to the general field of digital signal processing. 

[0004] More particularly, the present invention relates to multiresolution analysis for signals or data 

of any dimensionality using a non-separable, radial frame, multi-resolution analysis in 

multidimensions. 

2. Description of the Related Art 

[0005] Digital signal processing, and, in general, the manipulation of information of all types by 
digitization, is at the heart of the computer-based approach to a vast range of problems in science, 
engineering, technology, modern economics modeling, data storage and retrieval, etc. There exist 
many robust approaches for problems which are intrinsically one-dimensional, and the theory of 
how one systematically parses the information content into small, manageable "chunks" is well 
developed. The essential idea is that information can be characterized either in a "physical" or 
"measurement" domain which, for example, we take to be the "time domain", or in a 
complementary, mathematical domain referred to as the "Fourier domain" (which we shall refer to 
as the frequency domain) .It is most useful if there are natural limitations on how much information 
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in the measurement or time domain is required to characterize a given amount of information in the 
frequency domain. 

[0006] In one-dimension (1-D), the best possible situation is when, e.g., only a finite range or 
"band" of data in the frequency domain is needed to characterize completely the underlying 
mathematical behavior. Such a situation is said to be "band-limited" and if one can capture, without 
any loss, all the frequency components contained in the band, the signal or phenomenon is exactly 
captured. Further, if the signal becomes "contaminated" by extraneous signals with frequencies 
outside the band range, then these can be eliminated by sending the signal through a "filter" that 
strains out everything except frequencies in the physical band range. This is accomplished 
mathematically simply by multiplying the signal (plus noise) by a function that is 1 for frequencies 
within the band-limit and 0 for all other frequencies. Such a filter is fundamental in all areas of 
signal processing and information analysis, and it is called "ideal filter" of "ideal window". 
[0007] The amazing fact is that when one must treat signals that are of higher dimension including 
the simplest case of (two dimensional (2-D) data), the only rigorous way that exists to create such 
an ideal window (and to shift or translate them to different frequency bands can be captured leading 
to a multiresolution) is by multiplying 1 -D ideal windows for each degree of freedom. Such products 
are said to be "separable" and they are inefficient for studying data sets or signals for which there 
is not a natural separation of information content along orthogonal directions. It is expected that in 
the absence of knowledge of such a directional bias in the signal, the best approach would treat the 
data in the most isotropic manner possible. 

[0008] Thus, there is a need in the art of signal processing for the construction of improved multi- 
resolution analysis techniques for extracting information from complex scientific signals, especially 
processing techniques that involve non-separable, radial frame, multi-resolution analysis in one or 
more dimensions. 

SUMMARY OF THE INVENTION 
[0009] The present invention provides a signal processing technique using a completely isotropic 
ideal filter, and then translating and dilating the filter to create a multiresolution analysis. The 
technique includes the following constructed components: 1) isotropic ideal windows or filters in 
any number of dimensions; 2) translations and dilations to create completely isotropic low pass 
filters, high pass filters and/or filters that cover all frequencies and/or frequency ranges important 
to any appropriate signal processing problem; 3) isotropic scaling functions with translations of the 
low pass scaling function; and 4) associated wavelets that can be used to resolve a multidimensional 
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signal into various resolution levels, where the technique admits as many levels of multiresolution 
as desired, so that the high pass components can be further parsed into their own "relative low pass" 
and "relative high pass" components. 

[0001 0] The present invention provides scaling functions, wavelets and various isotropic non- 
separable ideal windows or filters that are the components needed to construct completely isotropic, 
intrinsically non-separable low pass and high pass filters, scaling functions, and wavelets that also 
have the optimum localization of information in a time-frequency description. 
[0011] The present invention provides an improved method to analyze signals, data, information, 
images as well as other items of any number of dimensions, both in terms of efficiency and 
complexity of problems that can be addressed. The present technique can be used to analysis: (a) 
data compression and storage for streaming video, seismic imaging, digital medical imaging of all 
types, etc.; (b) image and signal enhancement, denoising and analysis for medical imaging, seismic 
imaging, satellite imaging and surveillance, target acquisition, radar, sonar, pattern recognition and 
analysis, etc.; (c) volume rendering and segmentation, motion analysis, etc.; and (d) as a basis for 
digital algorithms for solving ordinary and partial differential equations in science, engineering, 
economics, and other disciplines. 

[0012] The present invention further relates to a computer or computer readable storage medium 
having implemented thereon software corresponding to the signal processing techniques of this 
invention. 

DESCRIPTION OF THE DRAWINGS 
[0013] The invention can be better understood with reference to the following detailed description 
together with the appended illustrative drawings in which like elements are numbered the same: 
[0014] Figure 1 depicts a plot of each of the sets of D + k r ; 

[0015] Figure 2 depicts a partition of Q into 29 sets due to intersecting Q with B{e 0 € i9 €^€ 3 ); 
[0016] Figure 3-6 depict sets of C, for / equal to 1, 2, 3, and 4, respectively; 
[0017] Figure 7(a) depicts that each C, overlaps Z>, 

[001 8] Figure 7(b) depicts all point in 7* at which the vector of the entries of the main diagonal; and 
[0019] Figures 8(a) and (b) depict sets C, and C 2 . 

DETAILED DESCRIPTION OF THE INVENTION 
[0020] The inventors have found that a new class of frame multiresolution analysis can be used to 
construct ideal filters. There are no other such totally isotropic, multidimensional, and non- 
separable scaling and wavelet functions and low and high pass filters currently in existence. The 
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present invention is unique. The first tests of the invention are directed to two dimensional (2D) 
images and data sets. 

[0021] If one could pass a signal through a set of ideal, isotropic windows or filters, each window 
passing a different band of frequencies, this would decompose the signal into subsets of information 
- we say the signal is resolved into a number of non-overlapping subsets and this is the basic idea 
behind "multiresolution analysis" of signals or digital information. 

[0022] We construct the new frame multiresolution analysis in any desired number of dimensions, 
whose core subspace is defined by the integer translations of a finite set of frame scaling functions, 
which have minimal support in the frequency domain. (The Fourier transform of the time domain 
scaling function is the characteristic function on a specific frequency range.) 
[0023] The most apparent representative classes of prototypes filters are a Frame multi-resolutions 
defined by a single scaling frame function, whose Fourier transform is a characteristic or "indicator" 
function of a hyper sphere (i.e., a disk in two dimensions (2D), a sphere in three dimensions (3D), 
hypersphere in higher dimensions (nD)), centered at the origin of the coordinate system having 
radius 1/2 (so that 0<R<l/2). Since this scaling function has this form in the frequency domain, its 
associated low pass filter in the frequency domain will also be such a characteristic hypersphere (or 
"ball") D/2, with radius 1/4. A fundamental property of the obtained wavelets and filters is that they 
are the first ever, truly intrinsically non-separable wavelets and filters that can be realized in any 
number of dimensions. 

[0024] Additionally, the Fourier transforms of the associated frame wavelets will be equal to the 
characteristic function of a "hyper-annulus" (a difference of two disks D/2 and D in 2D, a difference 
of two spheres D/2 and D in 3D, a difference of two hyperspheres D/2 and D in nD), modulated by 
a phase factor. 

[0025] For the 2D case, we create up to four frame wavelets; each of these wavelets defines a high- 
pass filter, or a similar wavelet like form. In the time (or complementary) domain, both the scaling 
functions and the wavelets can be constructed as linear combinations of Bessel functions. These 
constructions are currently based on dyadic dilations, but other dilations are equally usable only 
affecting a shape of the low and high pass filters. For example, dilations induced by quincunx 
matrices yield low and high pass filters of the same structure. 

[0026] A unique merit of our construction is that it leads to the definition and realization of low and 
high pass filters, which are intrinsically non-separable and isotropic to a degree that has never been 
achieved in the past and, which are, in fact, optimal. Our construction can be explicitly realized in 
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any finite number of dimensions. The resulting scaling functions are interpolating in all cases. 
[0027] The present invention can be implemented on any processing unit that is capable of 
executing instructions corresponding to the mathematical constructs and filters set forth in this 
application. Suitable processing units include, without limitation, analog processing units, digital 
processing units or mixtures or combinations thereof. These processing units are generally 
components of a computer of other device including a processing unit and peripherals capable of 
human interaction (keyboards and the like). Suitable computers include those manufacture and sold 
through out the industry based on chips from companies like Intel, Motorola, IBM, HP, Sun 
Microsystems, Cirex, AMD, or others and sold in computers manufactured and/or sold by companys 
such as Dell, Apple, IBM, HP, Toshiba, Sony, or similar computers. The processing units and 
computers incorporating them are designed to execute software under the control of an operating 
system. Suitable operating systems include, without limitation, the WINDOWS operating systems 
from Microsoft, the OS operating systems from Apple, the LINUX operating systems available 
from a variety of vendors, or other windowing operating systems. The techniques set forth in this 
application can be converted to software code in an number of convenient computer languages such 
as Fortran, C, C+, C++, or the like or newer programing languages geared to programming 
mathematical expressions. 
INTRODUCTION AND PRELIMINARIES 

[0028] Let //be a complex Hilbert space. A unitary system tf\s a set of unitary operators acting on 
H which contains the identity operator /on H. Now, let D be a (dyadic) Dilation operator 

(DJ)(t) = 2 n/2 A2t), feL\W) (1) 
and T k be a Translation operator defined by 

(TJ)(t) ^ - k\ fe L 2 (E"), k e Z" (2) 
We refer to the unitary system U DZ „ := {Z)T k :j e Z, k e Z n } as an n-dimensional separable Affine 

system. 

[0029] This system has been extensively used in wavelet analysis for the construction of separable 
wavelet bases. In fact, only a few non-separable wavelet bases have been constructed and all these 
examples were exclusively given in two dimensions. However, an important drawback of these 
families of wavelets is the absence of enough symmetry and differentiability and the absolute lack 
of isotropy. These examples were also given with respect to a small class of dilation operators and 
all of them are compactly supported in the time domain. Apparently the whole issue of designing 
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wavelet bases in multidimensions still remains a mostly unexplored area, full of challenges and 
revealing interesting and surprising results. 

[0030] The motivation for the present invention stems from the following elementary observation: 
The low pass filter corresponding to the scaling function of the Shannon MRA is the indicator 
function of the interval [-54, Vi). This function is even and minimally supported in the frequency, i.e. , 

its Fourier transform is of the form & - Xa > where A is a measurable subset of R. Keeping in mind 

that even functions are also radial (a function is radial if it depends only on the radial variable) one 
might wonder, what is the multidimensional analogue of even, minimally supported in the frequency 
scaling functions. 

[0031] This particular problem motivated us to introduce the radial frame multiresolution analysis. 
Our construction is based on a very general multiresolution scheme of abstract Hilbert spaces 
developed by Papadakis in [ 1 6] , namely the Generalized Frame Multiresolution Analysis (GFMRA). 
The main characteristic of GFMRAs is that they can be generated by redundant sets of frame scaling 
functions. In fact, GFMRAs encompass all classical MRAs in one and multidimensions as well as 
the FMRAs of Benedetto and Li ([6]). 

[0032] In this invention, we construct non-separable Shannon-like GFMRAs of L 2 (W) whose 
scaling functions are radial and are defined with respect to certain unitary systems, which we will 
later introduce. We also derive certain of their associated frame multiwavelet sets. Our construction 
is the first of its kind. Scaling functions that are radial have not been constructed in the past. 
However, certain classes of non separable scaling functions in two dimensions, with some continuity 
properties with respect to dyadic dilations or dilations induced by the Quincunx matrix only have 
been constructed in the past (e.g., [7, 12, 1 1, 8], [4]). All of them have no axial symmetries and are 
not smooth, except those constructed in [5], which can be made arbitrarily smooth, but are highly 
asymmetric. Another construction in the spirit of digital filter design, but not directly related to 
wavelets can be found in [1] and [18]. The latter construction and this of curvelets (e.g., see [19]) 
share two properties of our Radial GFMRAs: the separability of the designed filters with respect to 
polar coordinates and the redundancy of the induced representations. However, our construction in 
contrast to those due to Simoncelli et. al. and to Starck et al. are in the spirit of classical 
multiresolution analysis and can be carried out to any number of dimensions and with respect to a 
variety of dilation matrices. 

[0033] The merit of non separable wavelets and scaling functions is that the resulting processing 
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of images is more compatible with that of human or mammalian vision, because mammals do not 
process images vertically and horizontally as separable filter banks resulting from separable 
multiresolution analyses do ([20]). As Marr suggests in his book [13], the human visual system 
critically depends on edge detection. In order to model this detection, Marr and Hildreth used the 
Laplacian operator, which is a "lowest order isotropic operator" ([14]), because our visual system 
is orientation insensitive to edge detection. Thus, the most desirable property in filter design for 
image processing is the isotropy of the filter. Radial scaling functions for multiresolutions based on 
frames are the best (and, according to proposition 5, the only) type of image processing filters that 
meet the isotropy requirement. 
Definitions 

[0034] Before we proceed, we need a few definitions and results from M. Papadakis, "Generalized 
Frame Multiresolution Analysis of Abstract Hilbert Spaces, 2001" 

[0035] The family : i e I) is a frame for the Hilbert space H y if there exist constants A> B > 0 
such that for every x E H, we have 

ah 2 i z K*,*,)! 2 < *r 

16/ 

We refer to the positive constants A 9 B as frame bounds. Apparently for every frame, its bounds are 
not uniquely defined. We refer to the frame as a tight frame if A = B and as a Parseval frame if A 
= B = 1 . A frame {*,. : i e I) of H is called exact if each one of its proper subsets is not a frame for 
H. Riesz bases are exact frames and vice- versa. The operator S defined by 

S x ={(x,x t )} ieI xeH 
is called the Analysis operator corresponding to the frame {jc,. : i 6 I}. Using this operator, we can 
construct a dual frame {*',. : i 6 1} of {jc,. : / 6 1} by setting x\ = (5*5)'^. Then, for every jc e Hv/e 
have 

[0036] We are interested in unitary systems #of the form #= U 0 G, where U 0 = {U :j eZ) and G 
is an abelian unitary group. We will often refer to G as a translation group. Unitary systems of this 
form generalize the affine system. 
Definition 1 
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[0037] A sequence { Vj} )n of closed subspaces of an abstract Hilbert space H is a Generalized Frame 
Multiresolution Analysis of H if it is increasing, i.e., ^ c ^ for every j e Z and satisfies the 
following properties: 

(a) Vj=V(V 0 ), jeZ 

(c) There exists a countable subset B of K 0 such that the set G(B) = {g(f> :ge G,fie B} 
is a frame of V 0 . 

Every such set 5 is called a frame multiscaling set for {Vj) r Every subset C of K ; such that G(Q 
= {g^ : g e C, ^ e C} is a frame of W 0 = V t C\ V 2 0 and is called a semi-orthogonal frame 
multiwavelet vector set associated with {Vj}j- 

[0038] G(B 0 is the canonical dual of G(B) y where B ' = {(S*S)' 1 ^:^efi}, where S is the Analysis 
operator corresponding to the frame G{B). Likewise the canonical dual of G(Q is the family G(C% 
where C = {(S*S)'* y/e C). We refer to B ' and as a dual frame scaling set corresponding to B 
and to C as a dual frame wavelet set corresponding to C. 

[0039] If J? is a singleton, we refer to its unique element as a frame scaling vector and, if H = Z^R"), 
we refer to its unique element as a frame scaling function. We also let W } = U(W Q ), for every j e Z. 
Thus, if C is a semi-orthogonal frame multiwavelet vector set associated with the GFMRA {V } } p 
then the set {D>gy/: j e Z, g E G, y/E C) is a frame for H with the same frame bounds as the frame 
G(Q. 

[0040] In order to accomplish the construction of the frame multiwavelet sets associated with a 
GFMRA {Vj} p we need the following additional hypotheses. 
[0041] There exists a mapping a \ G -* G satisfying 

gD = Do{g\ for every g e G 

This particular assumption implies that a is an injective homomorphism and o(G) is a 

subgroup of <7. (See [9] for proofs) 

\G : o(G)\ = n < where \G : a{G)\ is the index of the subgroup o{G). 
[0042] As mentioned before in this invention, we will exclusively use multidimensional affine 
unitary systems. Before proceeding further with the construction, we need the following definition: 
Definition 2 
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[0043] An /i x n invertible matrix A is expanding if all its entries are real and all its eigenvalues 
have modulus greater than 1 . A Dilation matrix is an expanding matrix that leaves Z" invariant, i.e., 
A(Z n ) c Z". 

[0044] The previous definition readily yields the following observations: (a) all the entries of a 
dilation matrix are integers, because such a matrix leaves Z n invariant; and (b) the previous 
observation implies that determinant of A (detA) is an integer. 

[0045] The multidimensional affine unitary systems we are interested in are the systems of the form 
%G, where % is a cyclic torsion free group generated by a dilation operator D defined by 



where A is a dilation matrix and G = {T k : k e Z n ) and G is isomorphic with Z". Using the definitions 
of translations and dilations, we can easily verify TJ) = DT Ak , thus o{T^) = T Aky for every k e Z". 
Therefore, <ris legitimally defined, because A(Z n ) q Z". Apparently, the quotient group Glo(G) is 
homeomorphically isomorphic with Z n /A(Z n ). Thus we have \G : o(G)\ = \ Z n : A(Z")\ = \detA\. Now, 
we set q 0 = 0,p = \detA\ and we fix q r e Z", for r = 1, 2, ... ,/> - 1 so that 



[0046] The translation group G is induced by the lattice Z". Although our results will be obtained 
with respect to this particular lattice only, our methods can be easily extended to all regular lattices, 
i.e., lattices of the form C(Z W ), where C is an n x n invertible matrix. Using a traditional approach 
of Harmonic and Fourier analysis, we give the definition of the Fourier transform on L^W): 



We also reserve !7to denote the Fourier transform on L\R n ). In addition, we adopt the notation T" 
= [- ! /2, Vif and ^ = (i 6 / c N) for the sequence defined by 

fl if /= i 



[0047] Before proceeding, we need to include some final remarks on our notation. If A is a subset 
of a topological vector space, then [A] denotes its linear span and A' denotes the closure of A. 
Moreover, if B is a matrix (even an infinite one), then [B\ denotes the i-th column of B. We 
conclude this section with the characterization of the autocorrelation function of a set of frame 



Z)/(r) = |det4| 1/2 y(^), 



f€L 2 (R tt ) 



/«)= l Rn f(t)e^dt 9 £eR 



n 
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generators of a shift-invariant subspace of £ 2 (R fl ), of a set of functions {<f> t : / e I) such that { T k fy 
: / e I) is a frame for its closed linear span. 
Lemma 1.1 

[0048] Let / c N and {(f> k : k 6 1} be a subset of L 2 (l"). Define 

*/.*(<?) = I A(£+mM(f+m) kJel^eT" 

meZ n 

anda,(# = (a u (£, a w (0, ...). 

[0049] Now assume that for every £ e / the function £ - fa* (£)|| /2 is in L\V) and that the linear 

operators <B(£) defined for a.e?. £ e T" on [S k : k e 7] by the equation ®(§S k = a*(£) satisfy the 
following properties: 

(1) <D belongs to LlfT", £tf{I))\ i.e., <D is weakly measurable and for a.e. £ e T" the 
operator <D(£) belongs to and ||0|| w = essup{||<D(£)|| : £e T"} < «> 

(2) Let P(§ be the range projection of G>(£) a.e. There exists JS > 0 such that for every 
x E P(^)(i\I)) we have B\\x\\ z ||0>(£x||. 

Then {T k fy : / e I, k € Z"} is a frame for its closed linear span with frame bounds B and 

[0050] Conversely, if {7^, : / e /, * e Z"} is a frame for its closed linear span with frame constants 

B, C, then there exists O e I W (T", such that ||<D|| W ^ C also satisfying 

®{§tM = a M Me/, a.e. V 
and property (2). Finally, {T k ^ : / e 7, # e Z"} is a Parseval frame for its closed linear span if and 
only if <D(£) is for a.e. % an orthogonal projection. 

[0051] The function <D is also known as the Grammian of the set : / e 7} . 
RADIAL FMRAs 

[0052] In the this section, we will develop the theory of singly generated GFMRAs of L 2 (W) 
defined by radial frame scaling functions. We refer to these GFMRAs as Radial FMRAs. According 
to lemma 1.1, the Fourier transforms of frame scaling functions cannot be continuous. Thus, such 
scaling functions cannot have a variety of forms, but this drawback can be rectified by using frame 
multiscaling functions. However, in this invention, we will be exclusively using Minimally 
Supported in the Frequency (MSF) frame scaling functions. A function is MSF if the modulus of 
its Fourier transform is of the form Xa* where A is a measurable subset of 1". 
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[0053] Our t ranslation group i s g roup-isomorphic toZ w ,sowecane asily s ee t hat t he r egular 
representation of G on t 2 (G), defines a group, which we denote by G* and is homeomorphically 
isomorphic to the discrete group Z n . Therefore, the dual group G* is homeomorphically isomorphic 



Since 0(£) = for every £ e T", we have that {T k $ : k e Z"} is a Parseval frame for its closed 
linear span, which from now on, we will denote with V 0 . We will consider dilations induced by 
dilation matrices A satisfying the following property. 
Property D 

[0055] There exists c > 1 such that for every jc e R n we have c\\x\\ z \\Ax\\. 
[0056] Property D implies \\A' { \\ * c 1 < 1. However, it is interesting to note that Property D cannot 
be derived from the definition of dilation matrices. This fact can be demonstrated by the following 
example. Let 



We can clearly see that A is invertible and leaves the integer lattice invariant, because all its entries 
are integers. However, 



Since one of the entries of A A has absolute value greater than 1, we get > 1, so A does not 
satisfy Property D. 

[0057] Now, define ^=Zy(F 0 ), wherey e Z. We will now establish K/C K 0 . First, let B= A T , where 
the superscript T denotes the transpose operation. Since (A 1 )' 1 = {A' 1 )* and the operator norm of a 
matrix is equal to the operator norm of its transpose, we obtain that dilation matrices A satisfying 
Property D, therefore, satisfies \\R X \\ < 1. Thus, iT^D) is contained in D. Next, we let be the 

measurable function defined on W such that = X B -\ ( rxX€) > for ever y £ e T "> which is 



to the n-dimensional torus T". So, instead of using G , we use T, recalling that we identified T" 
with the product space [-J4, Vi)". 

[0054] Now, let D be the sphere with radius Yi centered at the origin, and ^ be such that (f) = % D . 



A = 



( 2 5> 
.0 2) 
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periodically extended on l n with respect to the tiling of l n induced by the integer translates of T". 
Then, ^ belongs to L\T n ) and satisfies 



turn establishes V A c V Q and thus Vj c V^, u for every integer / Since ^Vj) = Z, 2 (2P(D)), for all j e Z, 

we finally obtain that both properties of the definition of a GFMRA are satisfied. From the 
preceding argument, we conclude that {Vfo is a GFMRA of L\W), singly generated by the radial 
scaling function ft So {Vj}j is a Radial FMRA of L\R n ). We may also occasionally refer to ^ as a 
Parseval frame scaling function in order to indicate that {T k 0 : k e Z"} is a Parseval frame for V Q . 
[0058] Following previously developed terminology and the notation, the analysis operator S 
induced by the frame scaling set maps V Q into L\T) and is defined by 



where e k (& = e lni ^ k) for every £e M". Since $ is a Parseval frame scaling function, we obtain that 
S is an isometry. Moreover, it is not hard to verify that the range of S is the subspace L 2 (D). 
[0059] Thus, the low pass filter m 0 corresponding to ^ is given by m 0 = SD Because we consider 
{ Vj}j as singly generated, we have only one low pass filter, .so M 0 , the low pass filter associated with 
the frame multiscaling set is equal to m Q . Since S is an isometry, we obtain Y = S, where Y is 
defined by the polar decomposition of S, namely S= F|£|.Infact, we have Y=S= 3\^. This implies 

that m Q = YD*0=m o , and, therefore, M Q - m 0 - M 0 . Let us now find m 0 . Taking the Fourier 
transforms of both sides of 



because <f){B%) - X B -\ (m (4) , for every £ e W. This implies that 2>V » belongs to V 0 , which in 




keZ 




(3) 



keZ n 



we obtain 




a.e. 



(4) 



Page 12 



Specification 

Express Mail No.: EL 990 553 971 US 



Robert w. strozier. p.l.l.c. 



We recall 

kBQ=MQk& a.e. (5) 

[0060] Unfortunately, the fact that, the set of the integer translates of # is not a basis for K 0 , but an 
over complete frame, does not automatically imply that \detA\ A ^ = /w 0 . However, both m Q and 
vanish outside D, so equations (4) and (5) imply 

m 0 ($=\detA\ l ' 2 z B - i{D) (&, £er (6) 

[0061] All radial functions of the form x& where D is a sphere centered at the origin with radius r 
< X A are radial Parseval frame scaling functions. We will not distinguish this particular case from the 
case r = 54, because the latter case is generic and also optimizes the frequency spectrum subject to 
subband filtering, induced by this particular selection of the scaling function 0. This frequency 
spectrum is equal to the support of the autocorrelation function of <j>, because every signal in V Q will 
be encoded by the Analysis operator with an C 2 (Z)-sequence, whose Fourier transform has support 
contained in D. Therefore, the frequency spectrum subject to subband filtering induced by {Vj}j 
equals D. This suggests that a pre-filtering step transforming a random digital signal into another 
signal whose frequency spectrum is contained in D is necessary prior to the application of the 
decomposition algorithm induced by { Vj}j. This pre-filtering step is called initialization of the input 
signal. In the light of these remarks, one might wonder whether we may be able to increase the 
frequency spectrum that these FMRAs can filter by allowing r > 1/2. We will later show that the 
selection r = 1/2 is optimum. 

[0062] The frame scaling function can be determined in terms of Bessel functions, because it is a 
radial function. 

m--^f, *>o (7) 

The proof of equation (7) can be found in [17] Lemma 2.5.1. 

[0063] We will not give any details regarding Bessel functions. However, the reader may refer to 
[ 1 7] and [3] for an extensive treatment of their main properties and of course to the bible of the topic 
[21]. Here, we only include the following formula. 
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The function, J a , given by the above equation is called the Bessel function of the first kind of order 
a. 

[0064] Apparently every function in V 0 is band limited, because its Fourier transform is supported 
on D. Because D is contained in T, we infer from the classical sampling theorem that if /is in V Q , 
then 

/ = I /(k)7> (8) 

keZ D V ' 

where the RHS of the previous equation converges in the L-norm and co(x h x 2 , > x n ) = ll — . 

q 

If P 0 is the projection onto K 0 , then applying P 0 on both sides of equation (8) gives 

/= Zf(K)P 0 (T k co)= £ /(k)r k />(*>) 

keZ B keZ n 

because P 0 commutes with the translation operator T k , for every k € Z". Since P 0 (&) = <j>, we 
conclude the following sampling theorem: 
Theorem 3 

[0065] Let/be in K 0 . Then, 

keZ n V ' 

where the RHS of equation (9) converges in the L 2 norm. Moreover, the same series converges 

uniformly to /, if we assume that /is continuous. 

Proof 

[0066] The first conclusion of the theorem has already been established. We will now prove the 
uniform convergence to / of the series in the RHS of equation (9) assuming that /is continuous. Let 
t e 1". Then, for N e N, we have 



/(0-I^/(W) 



(10) 
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dt 



01) 



(12) 



As N- °° the first term of the RHS of the previous inequality tends to zero. This establishes the final 
conclusion of Theorem 3. 



[0067] The continuity hypothesis that we imposed on /in order to derive the uniform convergence 
to / of the series in the RHS of equation (9) is not at all artificial. It is well known that since /is 
band-limited, / is almost everywhere equal to an infinitely differentiable function, namely the 

inverse Fourier transform of / . Thus, instead of using/itself, we can use the reflection of 
Remark 2 

[0068] Although <f> is a radial function, its dilations , for j * 0 may cease to be radial, for if j = 
- 1 , then 3(D*$) = \detA\ h X B -\^) m ^ ma y not ^ e 311 isotropic domain. However, in several 

interesting cases of dilation matrices A all the dilations of ^are radial. 
[0069] The preceding remark motivates the following definition: 
Definition 4 

[0070] An expansive matrix A is called radially expansive if A = aU, where a > 0 and Uis a unitary 
matrix. 

[0071] Expansive matrices satisfy a n = \detA\ and \\A\\ = a and apparently radially expansive dilation 
matrices satisfy Property D. When this is the case, we immediately obtain that all D 7 ^ are radial 
functions as well, and, in particular, 



Remark 1 




R> 0 



(13) 



Combining equations (4), (6) and (13) we conclude 
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Proposition 5 

[0072] Let A be a radially expansive dilation matrix, and D r be a sphere having a radius r centered 
at the origin. Then, there exists an r 0 > 0 such that, if r > r 0 and <f>- D'\ Xo r \ th en no measurable 

Z"-periodic function satisfies 

for a.e. £ in E". Thus, such a ^ cannot be a frame scaling function. 
Proof 

[0073] Let r > 0 and <f>= Xd,- Assume A = aU, where a > 1. Then, (f>(B%) = % D (£) a.e. £in 

r 1 

R n , which in conjunction with equation (14) implies M£) = , (£) for <a.e. £e If — > , 

then //(£) = 1 for a.e. £in T", which, due to the Z"-periodicity of//, implies = 1 ae. in R". This 

# a 
contradicts equation (14). Thus, r < -j=. Now, pick such an r. If r < r 0 = — — , then Z) r , and 

k + D rla , for every k e Z\{0} do not intersect. 



u 

[0074] Now, assume r > r 0 . In this case, we have 1 / 2 < r 0 < r < , N ex t translate T" by u = (1, 0, 

0, ... 0). Due to the periodicity of //, we have //(£)= 1 for a.e. £in the intersection of the sphere 

\r 11 

u + D , andu + T". Because r>r 0 , we can find jc such that maxi — — f < x< r n . Then, there 

exists a ball centered at (x, 0, 0, ... , 0), which is contained in the intersection of u + D r/a , u + T" and 
W\D r so, equation (14) fails to be true for every point in this ball. 

[0075] If V 2 < r < r 0 , then (f>- F 1 (Zd ) 

is a frame scaling function. This can be shown by 

invoking lemma 1.1, which establishes that {T k # : k E Z n } is a frame (but not a Parseval frame) for 
V 0 and the argument showing. V } = D'\L\BfD)). We will omit the details of this proof since we 
think that this particular case is not as interesting as the case r <, J4, because the FMRAs defined by 
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such frame scaling functions ^ still cannot filter the entire n-dimensional torus T". Having finished 
this intermezzo, we return to the initial hypothesis, r = X A. 

[0076] Let us now discuss the construction of certain frame multiwavelet sets associated with { Vj} r 
The cardinality of the frame multiwavelet sets associated with the same GFMRA may vary. This 
observation indicates that there is room for alternate constructions of GFMRA frame multiwavelet 
sets. However, all these sets must satisfy certain necessary and sufficient conditions, which we 
present in theorem 6. 

[0077] In the discussion that follows, we present two constructions of frame multiwavelet sets 
associated with { Vj}j. Each one has its own merit. The first one does not depend on the dimension 
of the underlying Euclidean space R", and we believe that it is the most elegant set of them all. The 
second one specifically applies only if the underlying space is I 2 and the dilation operators are 
defined by A = 2I 2 or 

'1 1 
- 1 1 



It can be seen that both matrices are radially expansive dilation matrices. The latter of these two 
matrices generates the so-called Quincunx subsampling lattice. Subsampling lattices are used in the 
applications of the Decomposition and Reconstruction algorithms. 
First Construction 

[0078] Using a theorem relating to the radial FMRA {Vj} p we first set Vj = F{Vj) and 
Wj = F (Wj ) , where jeZ. Recall that V 0 = F(V Q ) = L 2 ( D ) , and that the Fourier transform is 
a unitary operator on L*(BL*). Combining these facts with V_ x = L (B~ (D)) , we conclude 

W_ x = V o nv_\ = L 2 (Q) 

where Qis the annulus D H (2r ! (D)) c , and the superscript c denotes the set-theoretic complement. 

Since an arbitrary orthogonal projection R defined on a Hilbert space H maps every orthonormal 
basis of H onto a Parseval frame for R(H) ([2, 10]), we obtain that the orthogonal projection defined 
on L\V) by multiplication with the indicator function of Ogives a Parseval frame for L\ Q, namely 

the set {eyXo : keZ"}. 

[0079] Next, observe that each k e Z n belongs to exactly one of the elements of the quotient group 
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Z n /A(Z tt ); thus there exist a q and r 6 {0, 1, ... j? - 1} such that k = q, +A(q). Therefore, 
e k = e q ^A(q) • ^ e now define following functions: 

K = e qr%Q r 6 ((U p-l) 05) 

Apparently {e A(k) h r : k 6 Z n , r = 0, 1, . . . , p - 1 J is a Parseval frame for L\ Q, thus, for 

as well. Therefore, {t^F " 1 // r :k € Z V = 0,1,...,/?- l} is a Parseval frame for W mU because the 

Fourier transform is unitary. Setting y/ n =DT x h r (r = 0, 1, ... y p - 1), we finally have that {T k y/ r : 

k e V , r = 0, is a Parseval frame for therefore {y/ r : r = 0, is a Parseval 

frame multiwavelet set associated with the FMRA {Vj}j. This concludes the first construction of a 
frame multiwavelet set associated with {Vj)y 

[0080] The reader might wonder whether it is possible to give a more explicit formula for the frame 
wavelets y/ r . In the light of remark 2, % may not be radial as well. This may yield a rather 
unattractive time domain formula for all these wavelets. It is worth mentioning that y/ n where r > 
0, are never radial if % is radial. However, if A is a radially expansive dilation matrix and a = \\A\\ 9 
then 

IF - {2Ry n {2aRy n . R > 0 

Therefore, under this assumption, % is radial and 

a nl2 J nl2 {mR)- J.J71R) 

and forr = 1, 2, ... 

^ (0 = DT q D*y/ 0 (t) = ¥o (t- A~\) 

-'aJb 





A 




)- J„, 2 (4-^-'q r |) 


(2a 


/I 


r'qj)"' 2 
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It can be seen that in this case p = \detA\ = <f. 

[0081] We now continue with the preliminaries of the second construction. From now on and until 
the end of the present section, we work with GFMRAs of L 2 (R 2 ) only. 

[0082] One of the instrumental tools of this construction is the square root of the autocorrelation 
function <D, which is defined by A{£) 2 = a.e. on T. Also, the inverse ofA(g) is defined on the 
range of 0(£) and is denoted by A(£y x . It can also been proved that the range projection P of the 
Analysis operator S is defined by Pco(%) = P( where a) € L 2 (J n \ and that for a.e. £ e T 2 the 
range projection of 0(£) is the projection P(£). For the sake of completeness, it must be noted that 
P( ) is a projection- valued weakly measurable function defined on T 2 . Since 0 = % D , we deduce 

= Xo (£) a - e - m T 2 - The latter observation in conjunction with the preceding argument 
implies that A(& 1 = 1, if £e D. For all other £ 6 T 2 , we have A(§ = 0, so for these £ we adopt the 
notational convention A(gf x = 0. Last but not least, an abelian group very instrumental in the 
discussion that follows is the kernel of the homomorphism p defined by 

p(£)(k)= e ««* k >, keZ" 

The latter equation implies that, for every £ e T 2 , is the unique point in T 2 , such that + k 
= A T ^ . The kernel of p is homeomorphically isomorphic to a dual group of the quotient group 
Z 2 /A(Z 2 ) as shown in greater detail in reference [16]. Now, let us fix k r , where r = 0, 1, ... ,p - 1, in 
T 2 , so that Ker/?= {k, : r = 0, 1, ... ,p - 1}. 
Theorem 6 

[0083] Let / c N. Assume H : T 2 -> S(7 2 (/),C) . Define 



r=0 



Moreover, assume that the following conditions are satisfied 

(a) For a.e. e suppP 2 , where P 2 (£) is the range projection of the operator Q 2 (& , the 

operator Qi^A^eu)) :P 2(3(^(')) -> ^ 2 (^ 2 (/)) vertible and the 



functions 



&(^)l/> (0( ^ (/)) 



. <f-> 
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t 1 * I 



bounded. 

(b) For a.e. T 2 the closed linear span of the columns of the matrix 

M 0 (<f+k,) ff(#+k.) 



is equal to P(<^)(C P ) , where 



p-i 



= Z ® ^(# + M in T 2 and 



r=0 



( C ) 0= z;;>o(^+ k r r#(^+ k r ) 



[0084] If we define 



V, = I < n DT x "T 2 y 



where ja^ :z e 7,m, ft e z| are defined by the equation 

[#(•)], = I 

then, { \f/ i : i e I) is a frame multiwavelet set associated with the FMRA {^} ; . 

[0085] A measurable, Z 2 -periodic operator-valued function H , satisfying the hypotheses of the 
previous theorem is called a high pass filter associated with M 0 . If the dilation matrix satisfies 
Property D, then one choice for H following from equation (15) is 

H = ( e q 0 %Q> e qi ZQ >•••>%_, Xq) 

Let us first study the case where the dilation matrix A = 2f. In this case, it is well-known that p = 
4 and, that we can set k, = ( ! /2,0), k 2 = {Vi.Vi) and k 3 = (O//2) and recalling that addition in T 2 is 
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defined modulo the integer lattice Z 2 . 
[0086] Thus, we have 



0 



0 0 0 



0 
0 

0 



On the other hand, according to theorem 6, we must first determine the values of P before finding 
the high pass filter H . All the values of P are 4 x 4 diagonal matrices whose diagonal entries are 
either equal to 1 or 0. Therefore, the range of P is finite. 

[0087] So, we can find a partition of T 2 , say j B (€q ^ ^ ^ ) j where (5,, e l9 £,) is the vector formed 

by the entries of the main diagonal of an arbitrary value of P . Each s p , where p - 0, 1, 2, 3, takes 
only two values, namely 0 and 1 . 

[0088] Since each of the sets D + k r , where r = 0, 1, 2, 3, overlap with at least another one of these 
sets, there will be no values of P with a single non zero diagonal entry as shown in Figure 1. The 
definition of the addition operation on T 2 implies that D + k, is the union of the two half disks with 
radii 1/2 centered at k, and -k,; D + k 3 is the union of the two half disks with radii 1/2 centered at 
k 3 and -k 3 ; and, D + k 2 is the union of the four quarter disks with radii 1/2 centered at k 2 , -k 2 , 
(- 54,54) and { l A- X A). Since all four sets D + k r , where r = 0, 1, 2, 3, are symmetric with respect to 
both coordinate axes, it follows that all sets B^ ^ ^ share the same symmetry property. This 

observation contributes a great deal in identifying these sets. 

[0089] Referring now to Figure the subregions of T 2 corresponding to each one of the vectors (% 
£1, s^) are depicted. According to theorem 6, the values of the high pass filter H must be row 
matrices. This is justified by the fact that {Vj}j is generated by a single scaling function. However, 
the range of every P(£) is a subspace of C\ Furthermore, according to hypothesis (b) of theorem 
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**** A 

6, the columns of the modulation matrix must span P(£)(C ) . Thus, we anticipate that the 
modulation matrix must have at least three more columns. So, must be at least 1 * 3 matrix. 

For reasons that will become more clear herein, we choose H{£) to be 1 * 4 matrix, for every £e 
T 2 , namely 

The factor 2 in the RHS of the previous equation is a normalization factor that helps to obtain a 
simple form for the each of the functions h g . According to the conclusion of theorem 6, the columns 

of H , i.e., the functions h s (/ = 1, 2, 3, 4) , define a frame multiwavelet set associated with { 
[0090] Therefore, the modulation matrix has the following form: 

'WS £(0 Kit) y 

Zo, Ml i$ 5(£+k.) £(#+k|) . 2 

Z D /2 +k2 (0 M£+k 2 ) ^(£+k 2 ) ^(£+k 2 ) £ 4 (£+k 2 ) fl ^ mT 

Ud/m 3 (^) M£+k 3 ) %(£+k 3 ) ^(£+k 3 ) /T 4 (£+k 3 ) ; 



The disk D/2 has radius 1/4, so this disk and all its translations by k r (r = 1, 2, 3) have null 
intersections. Thus, for every £in T 2 , the first column of the modulation matrix has at most one non 
zero entry. Since, for every £in T 2 , the columns of the modulation matrix must span ?(£)(C 4 ) , we 
obtain the remaining columns of the modulation matrix, so that together with the first column they 
form the standard orthonormal basis of P(£)(C 4 ) . Thus, the high pass filters h t are the Z 2 -periodic 

extensions of the characteristic functions of certain measurable subsets of T 2 . Next, we will identify 
those subsets of T 2 , which we will denote by C i9 where i = 1 , 2, 3, 4. 
Remark 3 

[0091] Let Q be the first quadrant of T 2 . Then it is not difficult to verify that the family {Q + k f : r 
= 0, 1, 2, 3} forms a partition of T 2 , in the sense that T 2 = U* =0 Q + k r , but the intersections of 

every two of the sets Q + k r have zero measure. Now, let £be in Q + k„ then £= £ 0 + k n where £ 0 
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E Q. Without any loss of generality, we can assume r = 1 . Then, 

K)MZ+ k 2 )^.(# + k 3 )) r = + kOA(£U(^ + k,),*i(6 + k 2 )) r 

Thus, the values of the modulation matrix are completely determined by its values, when granges 
only throughout the first quadrant. 

[0092] As we have previously mentioned, the family £x €2 £ ^ j > where (6^, e u ranges 

throughout the vectors formed by the diagonal entries of the values of P , is a partition of T 2 . 
Therefore, jgf] ^( fo , M2 ,ff 3 )} is a partition of Q. Furthermore, each of the sets Qf] B(e 0t e lt e 2t e 3 ) 

is, in turn, partitioned into a finite number of subsets which are formed by the intersections of 
QC\^(s 0i £ U € 2> e^ with each one of the disks y, — + k„ y + k 2 , y + k 3 and the complement of their 

unions. This results in a partition of Q into 29 sets as shown in Figure 2. We denote these sets by 
E s9 where 1 ^ s <> 29. We now have to obtain, for every £in each of the sets E s , the remaining four 
columns of the modulation matrix, so that they span P(£)(C 4 ) . This process is not difficult to carry 

out. However, for the sake of the clarity, we deem it necessary to show how to specifically 
accomplish this task, when £ belongs to three of the sets E s . 
Case s = 1 

D 

[0093] This set is contained in the complement of the union of the disks — + k r , r = 0, 1, 2, 3, so 

the first column of the modulation matrix at £is equal to zero. Now, let ^e^. On the other hand, 
P(£)(C 4 )=C©0©C©C,sowe choose to complement the modulation matrix by setting its 

second, third and fourth columns equal to (1, 0, 0, 0) r , (0, 0, 1, 0) r and (0, 0, 0, l) r , respectively, and 
the fifth column equal to zero. 
Case s = 18 

[0094] Let £ e E n . Then, P(£)(C 4 ) = C©C©CeC. Moreover, <f = y . This suggests the 
following form for the modulation matrix at £ 
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( \ 0 0 0 0^ 

0 10 0 0 

0 0 10 0 

^0 0 0 1 0; 



Case s = 19 

~ . D 
[0095] Let£eE l9 . Then, P(£)(C 4 ) = 00 C0 C 0 O.Now, £e (—+ k 2 ) yielding the following 

form for the modulation matrix at £ 

r 0 0 0 0 0^ 

0 10 0 0 

1 0 0 0 0 
^0 0 0 0 0; 



It is now easy to verify that C„ where i = 1, 2, 3, 4, are the sets depicted in Figures 3, 4, 5 and 6, 
respectively. Let us now briefly review the case 

i \\ 



A = 



-1 V 



In this case,/> = |det4| = 2. It is also not hard to verify k, = (54,54) and 



Thus, D + k, is now the union of the four quarter disks with radii 54 centered at the vertices of the 
fundamental domain T 2 . 
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[0096] Each disk overlaps with D as shown in Figure 7(a). This overlapping, as in the case of A = 
2/ 2 , yields a partition of T 2 , namely the collection of subsets B( £Qy€x ) where, (% £,) is the vector 

of the entries of the main diagonal of an arbitrary value of P , and B(e Qy e x ) contains all points in 

T 2 at which the vector of the entries of the main diagonal of P is equal to e { ) as shown in Figure 
7(b). The low pass filter is now given by 

mo(£)= V2^(^),^gT 2 

as shown in equation (6)). This expression of the low pass filter follows from the form of the dilation 
matrix which is a composition of a rotation by n/4 matrix and >/2/ 2 . We can now take 



for every £ e T 2 . The modulation matrix now has a much simpler form, namely 

a.e. in T 2 



f xAZ) Kit) 

V2 



Zo (0 *.(£+k,) hit+K) 
v v2 +ki ; 



[0097] Let us now set Q to be the closed square whose vertices are the mid points of the sides of T 2 . 
It is not hard to see that Q + k, is the union of the four orthogonal isosceles triangles defined by the 
vertices of Q and T 2 . Obviously, {(?, Q + k,} is a partition of T 2 modulo null sets. An argument 

similar to the one in remark 3 shows that it is enough to determine the filters h k (/ = 1, 2) only 

on Q. It will also be helpful to observe that the sides of Q are tangent to the circle of radius 2 

centered at the origin and that Q can also be partitioned by the sets Qf) B {Sq €]) , where (% *i) = 

(1,0), (1,1) as shown in Figure 7(b). Each of these two sets will also be partitioned by its 
intersections with each one of + k, and the complement of the union of the latter pair of sets 

as shown in Figure 6. This, now results in a partition of Q into 17 sets. 
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[0098] Arguing as in the case of A = 2/ 2 , we can now obtain the sets C, and C 2 , so that 
K (£) = Zq (£) > where £ E T 2 and i = 1, 2 as shown Figure 8(a) and 8(b), respectively. 
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[00100] All references cited herein are incorporated herein by reference for all purposes 

allowed by law. While this invention has been described fully and completely, it should be 

understood that, within the scope of the appended claims, the invention may be practiced otherwise 

than as specifically described. Although the invention has been disclosed with reference to its 

preferred embodiments, from reading this description those of skill in the art may appreciate changes 

and modification that may be made which do not depart from the scope and spirit of the invention 

as described above and claimed hereafter. 
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